Heat Transport in Confined Strongly Coupled 2D Dust Clusters 
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Dusty plasmas are a model system for studying strong correlation. The dust grains' size of a few 
micro-meters and their characteristic oscillation frequency of a few hertz allows for an investigation 
of many particle effects on an atomic level. In this article, we model the heat transport through an 
axially confined 2D dust cluster from the center to the outside. The system behaves particularly 
interesting since heat is not only conducted within the dust component but also transfered to the 
neutral gas. Fitting the analytical solution to the obtained radial temperature profiles allows to 
determine the heat conductivity k. The heat conductivity is found to be constant over a wide range 
of coupling strengths even including the phase transition from solid to liquid here, as it was also 
found in extended systems by V. Nosenko et al. in 2008 pj. 



I. INTRODUCTION 



Strongly correlated systems are of high interest in mod- 
ern physics. Strong correlations can emerge at very dif- 
ferent physical regimes. In dusty plasmas, the high par- 
ticle charges are responsible for the strong correlations, 
while ions in traps or laser plasmas can be strongly cor- 
related due to the small particle distances 2J. 

Besides static properties, experiments and Molecu- 
lar Dynamics simulations provide insight to dynamic 
properties. Particle transport (diffusion) is well inves- 
tigated [3, 4 . Also, the heat transport is is an impor- 
tant property that was experimentally investigated in 
extended systems [I] [S] . One question is how the strong 
correlations affect the transport properties. 

In this article, we investigate the heat transport in a 
finite 2D cluster in an experiment-oriented simulation. 
After a brief presentation of the system model and the 
simulation technique in Sec.lnj we present simulation re- 
sults for different parameters in Sec. |III| The analysis 
focuses on the temperature of the dust particles' random 
motion and the radial temperature profiles under differ- 
ent heating conditions. These profiles are compared to 
the solution of the stationary heat transport equation for 
the 2D cluster in Sec. |III A[ which allows to derive the 
thermal conductivity k under certain assumptions. Fi- 
nally, we discuss the results for k and the dependence of 
k on heating power, equilibrium coupling strength and 
screening parameter in Sec. |IV| 



II. EXPERIMENTAL SETUP AND 
SIMULATION MODEL 



2D dust clusters are usually realized by using a lower 
electrode with a cavity. The plasma sheath which con- 
tains the dust grains then reflects the shape of the cav- 
ity. This shape allows for a horizontal confinement of the 
dust. In order to investigate heat transport, it is neces- 
sary to heat to the dust component in a restricted area. 
The laser manipulation technique [51 |7j is well suited for 
this purpose. Several randomly moving laser spots ac- 
celerate dust particles in different directions due to the 
momentum transfer by the radiation pressure. A spa- 
tially inhomogeneous heating is achieved by restriction 
of the area that is scanned by laser beams. 

The unheated A^-particle system is described by the 
Hamiltonian 
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is the pair distance and Ad is the 
Debye length. This model assumes that all particle are 
equal in mass m and charge Q. Wake effects that oc- 
cur in streaming plasmas are neglected, since a single 
layer of particles is subject to our studies. The horizon- 
tal interaction is adequately describe by a Yukawa poten- 
tial [F . This Hamiltonian is transformed to a dimension- 

less form by mtroducmg to = 1/w, Iq = (^j^^j^ 
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^0 = I 74 — Y^ I ^^ units for time, length and en- 
ergy. The Debye screening is then characterized by the 
screening parameter k — Iq/Xb- 

While the dust component is treated exactly, electrons, 
ions and neutral are treated statistically in the Langevin 
model. In this model, the amplitude of the stochastic 
force is expressed in therms of the friction coefficient 7 
between the dust particles and the neutral gas and the 



equilibrium temperature T^q. T^q and 7 are input pa- 
rameters to the simulation. We use Langevin Molecular 
Dynamics simulations [^ to propagate the equations of 
motion for the dust component and obtain the trajecto- 
ries. The heating lasers are included as space and time 
dependent forces [7] 
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in this model. The index / counts the laser. The force 
amplitude is given by Pq, a^/y describe the anisotropic 
spot profile and e; the a unit vector in beam direction. 
r/(i) = {xi{t),yt{t)) is the trajectory of the moving spot 
in the levitation plane. We use two pairs of laser spots, 
one pair in x and one in y direction. The lasers of each 
pair accelerate the particles in opposite direction, thus 
the average total momentum transfer vanishes. The laser 
trajectories are generated randomly. The spot I moves 
across the levitation plane with a velocity v;. When it 
reaches a border of the heated window, the motion is 
reverted in this direction and a new speed is randomly 
chosen within an interval. This method was shown to al- 
low for a controlled heating of the dust component with- 
out altering the properties of the plasma and the neutral 
gas [i[7]. 

We use the Coulomb coupling parameter F ~ 
Q"^ I {AneolokBT) [101 HI] as ration of typical interaction 
energy and thermal energy to characterize the coupling 
strength. As characteristic particle distance, the unit 
length l(j is used. 



III. SIMULATION RESULTS 

Fig.Illshows an example for the trajectories in a cluster 
with N = 200 particles, with k = 1 and Toq = 200. The 
four laser spots move within a square with a side length 
a = 1 in units of Iq. While the motion is restricted to a 
small region around the equilibrium position for particles 
at the cluster's edge, particles in the central region are 
less localized. 

The trajectories in Fig. [l] show that the inner parti- 
cles move greater distances within the same time than 
the outer ones. In order to quantify the associated tem- 
perature gradient, we calculate the radial temperature 
profile. Therefore, the simulation space is divided into 
circular rings of width Ar. The width has to be chosen 
small enough to obtain a sufficient radial resolution and 
large enough to obtain a sufficient number of data points 
per ring. Ar = 0.3 in units of the trap length turned 
out as a proper choice. The temperature for the ring is 



calculated as k^T 



((v^ 



(v) j by taking the av- 



erage over the whole simulation after an equilibriation 
phase of typically icq — 20w~^. The temperature profile 
reaches a steady state a few plasma periods lu~^ after 
powering the heating lasers in the simulation. The time 
resolved temperatures for three different radii are shown 



in Fig. [2] Since the dust cluster is trapped, the collective 
fiow velocity (v) = vanishes. Figure Is] shows the veloc- 
ity distributions of different rings from the center to the 
outside of the cluster. The distributions are Maxwellian 
for the inner rings as well as for the outer rings. The 
width of p{v) is large in the center and decreases mono- 
tonically towards the cluster's border. That means, the 
temperature decreases with the radius, as expected. 



A. Analytical model 

In this section, we present an analytical model that 
can describe the radial dependence of the temperature. 
Our goal is to describe the dust temperature outside the 
central heated region. Therefore, we consider a incom- 
pressible fluid model and add the thermal coupling to the 
background gas with temperature Teq in the heat trans- 
port equation |12j . 
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div (fcVT) + 5v - 27n(r - Teq) 
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On the right hand side of Eq. (|3| , we find the convective 
time derivative of the temperature. The right hand side 
consists of the thermal conduction, the shear heating as 
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FIG. 1. Trajectories in the first quadrant of a heated Yukawa 
cluster during At = lOto ■ Each one arbitrary particle is high- 
lighted in the central, middle and outer region. The gray 
pattern around the origin indicates the array scanned by the 
laser spots. 

Parameters: particle number A'' = 200, screening parameter 
K = 1, equilibrium Coulomb coupling parameter Fe, — 200, 
friction parameter 7 = 0.5 and heating power P = 90. 
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FIG. 2. Time evolution of the local temperature in the sim- 
ulation for different radii. The lasers were enabled at f = 0. 
While the temperature at the inside saturates on a time scale 
r = 1.7, the relaxation is slower for the outer particles with 
r = 4.3. The temperatures where average over 40 indepen- 
dent simulations. 
Parameters: A^ = 200, k = 1, Tcq = 200, 7 = 0.5, P = 50 




2 3 4 

radial coordinate r 

FIG. 3. Velocity distribution p{vx) of particles within con- 
centric rings (A'^ = 200, moderate heating power P = 40, 
equilibrium coupling strength V^q = 200). The data points 
are well fit by Gaussians (filled curves) with decreasing width 
CTi, = \Jk-BT{r)/m. The velocity distribution in j/-direction 
pivy) coincides with p(vx) and is omitted for this overview. 
The heated area is indicated by the gray pattern. 



function of the spatial derivative of the velocity vector 
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where the indices i and j take values 1 , 2 in two dimen- 
sion and the squaring implies summation over i and k. 
The last term describes the losses to the neutral gas. If 
we wanted to describe the central region where the laser 
heating takes place, we would have to introduce another 
source term. When we consider the dust cluster long 



dT 



partial time derivative ^ vanishes. 

The stationary heat transport equation[ 



reads 



Cpuyr ■ VTfcB = div(A:Vr) - 2-fn{T - Tcq)kB + S^ , (5) 

where Cp is the specific heat, n the spatial number den- 
sity and V is the flow velocity [T]. While the equilibrium 
temperature Tcq and the friction coefficient 7 are input 
parameters to the simulation, the heat conductivity k is 
unknown. We further neglect the convective heat trans- 
port on the left side and heat input by viscous conversion 
of energy form a shear flow ^v . This is motivated by the 
observation that no shear flow is observed in the simula- 
tion data. We check the validity of the assumptions later 
by comparing the temperature profiles in the simulations 
with the analytic solution. The reduced heat transport 
equation then reads 



div(fcVT) = 2jnkB{T - Toq) . 



(6) 



Due to the symmetry, we expect the temperature to de- 
pend on the radial coordinate r only, and rewrite Eq. m 
in polar coordinates as 
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When losses to the neutral gas are neglected (7 = 0) 
and the heat conductivity k is independent of r, the 
solution of Eq. ([7| is a logarithmic temperature profile 
T{r) — ciln(r/ro), with the integration constants rp 
and ci. As shown in Fig. |4j this solution does not re- 
produce the velocity profiles from the simulation. The 
temperature in the simulations decays much faster than 
the logarithmic profile. The logarithmic solution would 
decay indefinitely, while we expect the system to be al- 
ways warmer than the surrounding gas due to heating 
applied. Hence, losses to the neutral gas by friction are 
important for the spatial heat distribution. These losses 
are responsible for the faster decay. The general solution 
of Eq. ([7]) with a heat conductivity k independent of r is 
readily found and given by modified Bessel functions of 
first kind /g and second kind Kq, 
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where the abbreviation b = 2^nkB/k is used. Since the 
dust cluster has a finite radial extension R, no heat can 
be transferred further outside. A heat fiow equal zero at 
r = i? is ensured in the solution, when the temperature 
gradient vanishes. When Eq. (Is]) is differentiated, we find 
the second coeflicient B = AKi{VbR) / Ii{VbR) . 

The solution Eq. ([8]) then has two remaining free pa- 
rameters A and b which depend on the details of the 
laser heating, in particular, the heating power (see be- 
low). Since the heating is not included in Eq. ^ we 
determine these parameters by fitting the function to 
the simulation data. The central region where this input 
takes place is excluded from the fit. While the amplitude 
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FIG. 4. The simulation data (symbols in upper plot, A'^ = 
200, n — 1, Teq — 0.005, 7 = 0.5) are fit by both a logarithmic 
temperature profile T{r) = ci ln(r/ro) and by modified Bessel 
functions, Eq. ([8|. The latter fits the data well while the 
former solution does not. The central region was excluded 
from the fits since the power input by the lasers takes place 
in this region. The spatial density (lower plot) is alsmost 
constant in the central region where most of the heat loss to 
the neutral gas takes place. 



A is connected to the total power input, the thermal con- 
ductivity k is contained in b. The total heat that is trans- 
fered to the neutral gas outside the heated region (i?^) is 
calculated by integration of dPhcat — 2n'-fkB{T — Tcq)dA, 



Phcat = 27r / rdr 2n^kB{T{r) - T^q) 

jRt 
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Since the B/A — Ki{\/bR) / Ii{y/bR) decays fast, the sec- 
ond term can be neglected for clusters that are not to 
smaU. The ratio is B/A < 2 ■ IQ-^ for VbR > 5 while 
it is of the order of 1 at the inner radius with \/bR « 1. 
The total heating power is hence approxmimatly given 
by 



A ■ R / 
^heat ~ 47rn7 — yr^Ki i^VbR, 



(10) 



The dependence on the inner radius Rt remains, since 
Eq. ( 10 ) is the heat loss to the neutral outside the cen- 
tral region only and hence, depends on the size of this 
excluded region. 

Eq. ([9]) establishes the relation of the two free param- 
eters to the input power P. Since the radial temperature 
profile is well fit by the modified Bessel functions as so- 
lution of Eq. (l7|, see Figs. I4J [slandlBJ we conclude that 
Eq. (l7| correctly describes the physics. 



In this section, we investigate the infiuence of T^q and 
P on the thermal conductivity. A central question is, 
if the heat conduction changes during the transition be- 
tween the solid-like and liquid-like regimes. The melting 
point of a macroscopic 2D Coulomb (k = 0) system is 
well known as T = 137 \W . Due to finite size effects and 
finite screening parameter k, do not expect a sharp tran- 
sition temperature but a transition range |14| and choose 
the invested range of F accordingly. 

As a first parameter scan, the equilibrium temperature 
of the Langevin thermostat is varied at constant heating 
power. The radial temperature profiles for four different 
values of Tcq are shown in Fig. [5| We used a compar- 
atively small laser power for this plot, since the back- 
ground temperatures are of the same order of magnitude 
as the central temperatures in that case. The tempera- 
ture approaches Tcq towards the cluster boundary in all 
four simulations. Varying T^q has also an effect on the 
temperature of the heated particles in the center. The 
temperature is increased by an increase of Tcq. This ob- 
servation is not surprising, since heat input by collisions 
with the neutral gas (i.e. via the stochastic force) and 
heat input by the laser force add up. While the ampli- 
tude of the temperature increases with Teq the parameter 
b remains constant. However, the dispersion of the sin- 
gle measurements is large as seen in Fig. [7] In order to 
estimate the error of 6, we performed a series of 20 sim- 
ulation for each value Fgq = 50, . . . , 300. Each run uses 
the same parameters but different random numbers and 
initial particle positions. Mean value and standard devi- 
• ation of b are also plotted in Fig. |8] Within the accuracy, 
we find that b = 1.57±0.01 (P = 50) and b = 1.66±0.01 
(P — 100) is constant when Tcq is varied. A series of 20 
simulations was also performed with moderate coupling 
Fcq = 1 and P — 100. The density profile becomes wide 
at the cluster boundary but is only changed weakly at 
the inside where most of the heat is transfered and lost 
to the neutral gas. With that, the result for b = 1.5 ±0.7 
indicates the independence of thermal conductivity and 
equilibrium temperature over a wide range. 

In a second parameter scan, the laser power is varied 
while the equilibrium temperature is fixed corresponding 
to Fcq = 200. The temperature profiles for four exem- 
plary powers are shown in Fig. [6] The higher power input 
results in an increased temperature in the cluster inside. 
For small powers, the equilibrium temperature is reached 
at the cluster's outer radius. For higher powers, the tem- 
perature of the outer particles is about twice as high as 
the equilibrium temperature. Interestingly, the param- 
eter b shows a linear dependence on the heating power 
within the accuracy of the estimation, see right side of 
Fig.i 
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FIG. 5. Simulation data (symbols) and fits by modified Bessel 
functions (lines) for different equilibrium temperatures and 
constant heating power P = 20. The central region was al- 
ways excluded. On the right axis, the temperature is trans- 
lated into a local coupling strength. 



FIG. 7. Values of the parameter 6 for different equilibrium 
temperatures Fcq. Each symbol represents a single simula- 
tion. The error bars show the variance of the single estimates 
of b. 
Parameters: N = 200, k = 1, 7 = 0.5, P = 100 
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FIG. 6. Simulation data (points) and fits by modified Bessel 
functions (lines) for different heating powers. The central 
region was always excluded. The equilibrium temperatures 
corresponds to Feq = 200 in all cases. 



IV. DISCUSSION OF THE RESULTS 

In this section, we give an interpretation of the sim- 
ulation results and compare them vifith experimental re- 
sults for bulk systems. The parameter b as introduced in 
Eq. ([8]) has the dimension of an inverse length squared. 
Hence, L — 1/ yh can be interpreted as an characteris- 
tic length for the heat transport. For arguments greater 
than 2, the modified bessel function Ki is well approxi- 
mated by an exponential decay. The temperature differ- 
ence T{r) — Toq drops by a factor « 2.5 each L in radial 
direction. The mean value of 6 = 1.66 from the sim- 
ulations for different equilibrium temperatures. Fig. |8J 
translates into L = 0.78 which is approximate the average 
inter particle distance. Another quantity which is often 
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FIG. 8. Values of the parameter 6 for different equilibrium 
temperatures of the Langevin thermostat (left) and for differ- 
ent heating powers (right). Each symbol represents the aver- 
age from 20 independent simulations with the variance of this 
average as error bar. The horizontal dashed lines indicate the 
mean value of 6 = 1.57 ± 0.01 for P = 50 and b = 1.66 ± 0.01 
for P = 100 in dimensionless units. The power dependence is 
fit by a linear dependence b = bo-\-m- P with 60 = 1.42 ±0.01 
and m = 2.68 ± 0.18 x 10"^. 
Parameters: iV = 200, k = 1, 7 = 0.5 



calculated in experiments [TJ [HJ [16] is the thermal diffu- 
sivity X — k/(ncp)^ where Cp is the specific heat capacity. 
The N particle system has 2N spatial degrees of freedom, 
where the common rotation is not connected with poten- 
tial energy and 2N momentum degrees of freedom. The 
specific heat is approximated as cp = {AN — l)/2iV w 2. 
Using the value of cp = 2, we obtain a thermal diffusivity 



of 
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(11) 



The conversion from dimensionless units into SI units was 
done for melaniin particles with dp = 6/iTO diameter, a 
charge of Q = 10, OOOe, a trap frequency uj = 5.5s~^ and 
a friction frequency 7 = aj/2 as typical plasma parame- 
ters. This result is in reasonable agreement with exper- 
imental findings for x ~ 9mm^ /s and Xtr ~ SOrnrn^ /s 
(transversal), xi ~ ^Qmm? / s (longitudinal) for 2D com- 
plex plasmas [U [15] and x ~ 1mm/ s in a liquid 3D com- 



plex plasma [T5]. 

We have shown that the radial heat transport in a 
confined 2D Yukawa cluster is well described by a fluid 
model. In this model, we could find an analytical form 
of the radial temperature profile. Using this model, 
Langevin Molecular Dynamics results show that the heat 
conductivity is constant over the investigated tempera- 
ture range. This range includes the transition region be- 
tween liquid-like and solid-like cluster. 
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